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Abstract 

The Mohave Rattlesnake ( Crotalus scutulatus ) is a highly venomous pitviper inhabiting the arid interior deserts, grasslands, 
and savannas of western North America. Currently two subspecies are recognized: the Northern Mohave Rattlesnake (C. 
s. scutulatus) ranging from southern California to the southern Central Mexican Plateau, and the Huamantla Rattlesnake 
(C. s. salvini) from the region of Tlaxcala, Veracruz, and Puebla in south-central Mexico. Although recent studies have 
demonstrated extensive geographic variation in venom composition and cryptic genetic diversity in this species, no 
modern studies have focused on geographic variation in morphology. Here we analyzed a series of qualitative, meristic, 
and morphometric traits from 347 specimens of C. scutulatus and show that this species is phenotypically cohesive 
without discrete subgroups, and that morphology follows a continuous cline in primarily color pattern and meristic traits 
across the major axis of its expansive distribution. Interpreted in the context of previously published molecular evidence, 
our morphological analyses suggest that multiple episodes of isolation and secondary contact among metapopulations 
during the Pleistocene were sufficient to produce distinctive genetic populations, which have since experienced gene flow 
to produce clinal variation in phenotypes without discrete or diagnosable distinctions among these original populations. 
For taxonomic purposes, we recommend that C. scutulatus be retained as a single species, although it is possible that C. s. 
salvini, which is morphologically the most distinctive population, could represent a peripheral isolate in the initial stages 
of speciation. 
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Introduction 

Current genetic diversity in the western North American vertebrate fauna has been shaped largely by Miocene and 
Pliocene orogenesis and subsequent Pleistocene glacial cycles (Douglas et al. 2006; Riddle & Hafner 2006; Man- 
tooth et al. 2013; Myers et al. 2016). During glacial periods, geographic distributions of widely distributed taxa 
were fragmented and reduced, starting the process of divergence in isolation. However, during interglacial periods 
and after the last glacial maximum, populations have expanded into secondary contact, providing opportunities for 
gene flow among once isolated populations (Wood et al. 2008, 2013; Myers et al. 2016). Although these cyclical 
global events likely influenced evolutionary histories of co-distributed taxa in similar ways, the extent to which 
speciation has occurred across disparate lineages has been equivocal (Schield et al. 2015). Regardless, the complex 
but partially tractable histories and demographics of these groups makes them ideal candidates for studying the nu¬ 
ances of speciation in vertebrates. 

The Mohave Rattlesnake ( Crotalus scutulatus Kennicott 1861) is distributed continuously throughout the inte- 
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rior arid zone of southwestern North America from southern California southeast to southern Mexico. Two poorly 
defined subspecies are currently recognized: the Northern Mohave Rattlesnake (C. 5 . scutulatus Gloyd 1940) occurs 
from southeastern California to the southern Central Mexican Plateau, and the Huamantla Rattlesnake (C. s. salvini 
Gloyd 1940) occurs in the matorral region from the eastern passes of the Trans-Volcanic Axis of State of Mexico, 
Tlaxcala, and Veracruz, south into the semiarid upper Balsas Basin (Gloyd 1940; Klauber 1972; Campbell & Lamar 
2004). Over this expansive distribution, these snakes inhabit desert flatlands, desert scrublands, arid grasslands, 
mesquite savannas, and oak scrublands. The taxonomic history is convoluted because early researchers considered 
C. scutulatus to be an “annectanf ’ form, intermediate between C. viridis Rafmesque 1818 and C. atrox Baird and 
Girard 1853 (Klauber 1930; Gloyd 1940). Due to the extensive geographic variation in venom composition across 
populations of C. scutulatus, including populations with neurotoxic, hemorrhagic, and mixed venom phenotypes, 
the species has been subjected to considerable research focusing on the various roles that diet, climate, and neutral 
processes have played in shaping venom phenotypes (Strickland et al. 2018a,b; Zancolli et al. 2019). Collectively, 
these studies have implicated local adaptation, rather than neutral processes, as a dominant factor driving fme-scale 
spatial variation in venom composition across its range. 

Crotalus scutulatus also was the focus of recent molecular phylogenetic studies, which inferred four well- 
supported and divergent mitochondrial clades that correspond with populations in the Mohave-Sonoran Desert, 
Chihuahuan Desert, Central Mexican Plateau (hereafter Central Plateau), and southern Mexican matorral (currently 
recognized as C. 5 . salvini ), respectively (Schield et al. 2018, 2019; Strickland et al. 2018b; Fig. 1). Schield et al. 
(2018) further demonstrated that nuclear genetic data additionally support these four distinct lineages, and provide 
evidence for past and/or ongoing gene flow between all geographically adjacent lineages. Accordingly, despite ge¬ 
netic evidence for as many as four distinct lineages, the appropriate taxonomic status of these lineages remains an 
open question. 
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FIGURE 1. Mitochondrial phylogeny of Crotalus scutulatus lineages reproduced from Schield et al. (2018), with arrows 
showing inferred gene flow between adjacent clades. Grey bars represent hypotheses of species limits corresponding to (A) a 
four-species model where each mitochondrial clade is a species, (B) a two-species model where the deepest mitochondrial split 
coincides with species (i.e., “expanded salvini ” hypothesis), (C) a three-species model where the deepest nuclear splits coincides 
with species (i.e., an alternative “expanded salvini' hypothesis), (D) a two-species model where current subspecies are elevated 
to species, and (E) a one-species model where C. scutulatus is recognized as a single species. Different sizes of arrowheads 
indicate direction of asymmetric gene flow between adjacent clades. 


Flere we investigate morphological variation across the distribution of Crotalus scutulatus and evaluate our 
inferences in context of recent molecular studies (Schield et al. 2018; Strickland et al. 2018b) to evaluate the taxo¬ 
nomic status of the species. Because our purpose was to interpret phenotypic patterns within the context of previ¬ 
ous genetic studies, we considered five alternative species hypotheses (Fig. 1): (A) each of four putative genetic 
lineages represents a unique species; (B) a more conservative hypothesis whereby the deepest lineage divergence 
(at the boundary between Chihuahuan and Central Plateau populations) divides C. scutulatus into two species (i.e., 
an “expanded salvini ” concept); (C) an alternative “expanded salvini ” concept, but which further splits Mohave- 
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Sonoran and Chihuahuan lineages, in congruence with the deepest nuclear splits (Schield et al. 2018); (D) each 
subspecies (as currently defined geographically) represents a species; or (E) C. scutulatus is maintained as a single 
species. Based on previous analyses of genetic data, each of these hypotheses has some merit, but a more definitive 
taxonomic assessment would benefit from a perspective from morphological variation. 


Materials and methods 

Species delimitation criteria. Most species concepts assume some level of evolutionary independence between spe¬ 
cies lineages, which ultimately results from genomic divergence and accumulation of reproductive barriers (Coyne 
& Orr 2004; Wu & Ting 2004). Under the perception of species as meta-population lineages, prezygotic reproduc¬ 
tive barriers are not a strict requirement among recently diverged species, but are contingent properties that could 
signal divergence along with other properties (de Queiroz 1998, 2007). Therefore, some level of gene flow might 
occur between nascent species, but should not be sufficient to hinder their independence as separate evolutionary 
lineages (de Queiroz 2007; Feder et al. 2012; Nosil & Feder 2012). Phenotypic data are useful not only for aiding 
physical recognition of species, but they also coarsely surrogate for genetic divergence as long as phenotypes are 
heritable. Fiere we follow a multivariate phenotypic cohesion criterion for evaluating whether morphologically 
cohesive groups correspond geographically with genetic groups or combinations of groups (i.e., whether shifts in 
distributions of multiple variables correspond with phylogeographic breaks). Alternatively, phenotypic traits might 
show clinal or mosaic variation across clades or even discordant variation, where phenotypically distinctive groups 
do not correspond with phylogeographic patterns. The rationale for phenotypic cohesion of multiple characters 
follows the basic evolutionary model of Fisher (1918), where continuous or quasi-continuous traits of individual 
species should be normally distributed under conditions approximating Fiardy-Weinberg equilibrium (Templeton 
2006; Cadena et al. 2018). Therefore, mixed normal distributions within a sample would be indicative of multiple 
potential species. 

Phenotypic data. We collected a series of one qualitative, ten meristic, and eleven morphometric characters 
for 347 alcohol-preserved specimens of subadult and adult Crotalus scutulatus from throughout its distribution. All 
specimens examined are accessioned in natural history collections in the United States, Mexico, and Great Britain 
(Appendix 1) and, when available, museum abbreviations follow Sabaj (2016). We recorded latitude and longitude 
for specimens using VertNet, Arctos, and collections data whenever possible, and georeferenced specimens from 
qualitative locality data when coordinates were not available using Google Earth. We followed Klauber (1972) Gut- 
berlet and Flarvey (2002), and Campbell and Lamar (2004) for terminology of morphological characters, and pro¬ 
vide character descriptions and abbreviations in Appendix 2. Our choice of characters was based on previous studies 
of variation in the C. virids-scutulatus complex, as well as an attempt to minimize strong dependence between some 
standardized scale counts. We determined sex by scoring presence of hemipenes. Total length and tail length were 
measured to nearest millimeter using a string and a metric ruler; all other morphometric characters were measured 
with digital calipers to the nearest 0.1 mm. We obtained snout-vent length (SVL) by subtracting tail length (CL) 
from total length (TL), and then subtracted head length (F1L) from SVL to obtain body length (BL). 

Among rattlesnakes, Crotalus scutulatus tends to have comparatively low subdivision of head scales in both 
the prefrontal and parietal regions, and approaches species of Sistrurus in this respect (Gloyd 1940). Both C. s. 
scutulatus and C. s. salvini have few intersupraocular scales (usually a minimum of two spanning a row), and have 
a distinctive large crescentric scale at the medioposterior border of each of the supraoculars (Cardwell 2016; Figs. 
2A-B). With respect to subspecies distinctions, previous authors have indicated that C. s. salvini differs from C. s. 
scutulatus in having less subdivision of crown scales, distal tail bands that are similar in color to dorsal blotches 
versus distal tail bands that are notably darker, and dorsal blotches that lack a distinct border of pale white scales 
(Klauber 1930; Gloyd 1940; Campbell & Lamar 2004). For each specimen, we scored whether pale blotch borders 
were distinct, indistinct, or absent (Figs. 2C-E). 

Statistical analyses. Because we evaluated morphological variation in the context of previous molecular stud¬ 
ies, we first assigned specimens of Crotalus scutulatus to one of four putative genetic lineages defined by Schield 
et al. (2018) based on locality of origin (Fig. 3). For our full dataset, we summarized geographic variation in indi¬ 
vidual quantitative characters across each putative genetic lineage. Because we defined condition of blotch borders 
as a qualitative character, we analyzed this trait separately using a Pearson chi-squared test of independence, which 
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evaluated whether frequencies of individuals assigned to each category were independent across putative genetic 
lineages. For all multivariate statistical analyses, we used a subset of individuals for which we had complete mor¬ 
phometric and meristic data (n = 243). Prior to analyses, we log-transformed morphometric variables, and excluded 
SVL and TL in order to include F1L independently in analyses, and because these measurements were used to derive 
BL and CL. We first applied principal components analysis (PCA) using the correlation matrix to reduce dimension¬ 
ality of our dataset and to evaluate gradients of morphological variation across the range of C. scutulatus in the con¬ 
text of correlated shifts in character values. We retained for further analysis all principal components that explained 
more than 5% of the total variation in the dataset. We plotted PC factor scores for the first three PC axes to evaluate 
whether discrete breaks in morphology were evident among lineages. We also applied the same character set to 
linear discriminant analysis (LDA) to determine whether training a model on labeled data (i.e., putative genetic 
lineages) would provide better separation of groups in multidimensional space (i.e., ‘morphospace’). From LDA, 
we also extracted the posterior probability of assignment of individuals to each lineage and averaged these values 
as a summary metric of discrimination among groups (see Meik et al. 2018). Because ordination of morphospace 
is subjective regardless of methodology, we also used model-based clustering to determine the optimal number and 
assignment of groups present in our dataset without a priori assumptions of classification (Fraley & Raftery 2007). 
Model-based clustering uses Gaussian mixture models (GMM) and Bayesian Information Criterion (BIC) to evalu¬ 
ate models with different numbers of groups (K), each with different multivariate volumes, orientations, and shapes 
(e.g., ellipsoidal, spherical, diagonal). 



FIGURE 2. Oblique views of head scalation in (A) Crotalus scutulatus scutulatus (MVZ 275547) and (B) C. s. salvini (MVZ 
275545), and representative examples of dorsal blotch condition as scored in this study: (C) pale dorsal blotch borders distinct 
(MVZ 275562), (D) pale dorsal blotch borders indistinct (MVZ 102990), and (E) pale dorsal blotch borders absent (MVZ 
275545). 

To evaluate general patterns in phenotypic variation among lineages while simultaneously accounting for sex¬ 
ual dimorphism, we used a series of two-way ANOVAs with genetic lineage and sex as the independent variables 
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and PC factor scores for the first five PC axes as the dependent variable. We regressed factor scores for the first 
five PCs against latitude and longitude using least squares linear regression to determine whether correlated shifts 
in character distributions showed continuous variation or corresponded with geographic breaks between adjacent 
lineages. To more generally evaluate an “isolation-by-distance” model using our morphological data matrix, we 
further converted raw data into a Mahalonobis distance matrix and performed a mantel test against geodesic dis¬ 
tances using 99,999 permutations. Because relative head size is an important taxonomic trait in some rattlesnake 
groups (e.g., C. mitchellii complex; Meik et al. 2015) and also indicates potential shifts in diet, we further analyzed 
the scaling relationship between HL (log mm) and BL (log mm) across C. scutulatus lineages using ANCOVA after 
confirming homogeneity of slopes. We used the MASS package to perform LDA (Venables & Ripley 2002), the 
mclust 5.0 package (Scrucca et al. 2016) for model based clustering, StatMatch (D’Orazio 2015) and ade4 (Dray & 
Dufour 2007) for mantel tests, and all other analyses were performed using the base statistics module in R v.3.4.2 
(R Core Team 2017). 



FIGURE 3. Geographic distribution of Crotalus scutulatus, with localities of specimens examined indicated by green symbols 
(Mohave-Sonoran lineage), yellow symbols (Chihuahuan lineage), orange symbols (Central Plateau lineage), and black sym¬ 
bols (C. s. salvini lineage). Pie-charts indicate relative frequencies of specimens with pale blotch borders absent (black), pale 
blotch borders indistinct (grey), and pale blotch borders distinct (white) for each genetic lineage. 


Results 

Summary data demonstrated extensive overlap of ranges for most quantitative characters across genetic lineages 
of Crotalus scutulatus, with a weak trend for lower mean or modal values for meristic traits in C. s. salvini (Table 
1). Frequencies in condition of dorsal blotch borders varied significantly among lineages (X 2 = 43.602, df= 6, p < 
0.0001), with the highest frequency of distinct pale blotch borders in the Mohave-Sonoran lineage, and a gradual 
decrease in frequency toward the south, with specimens assigned to C. 5. salvini having 100% frequency of pale 
borders absent (Fig. 3). With respect to PCA, the first five axes collectively explained 74% of the total variance in 
our dataset, and were retained for further analysis. A scatterplot of the first three PC axes showed relatively strong 
phenotypic cohesion across all C. scutulatus samples, with visible differences in dispersion of specimens assigned 
to each lineage, but without clear distinction among lineages (Fig. 4A). Moreover, a scatterplot of the first three ca- 
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nonical axes from LDA demonstrated considerably more separation among specimens assigned to genetic lineages, 
but likewise did not separate lineages into discrete subgroups (Fig. 4B). Overall accuracy of LDA assigmnent of in¬ 
dividuals to genetic lineages was 88.9%, and average posterior probabilities of assignment of individuals was 0.90, 
0.85, 0.74, and 0.82 to Mohave-Sonoran, Chihuahuan, Central Plateau, and C. s. salvini lineages, respectively. Most 
misidentifications from LDA were from specimens assigned to geographically adjacent genetic lineages, and the 
relatively low posterior probability of assignment of individuals to the Central Plateau lineage reflected a moderate 
frequency of misidentification of specimens between the Central Plateau and Chihuahuan lineages. As a final as¬ 
sessment of phenotypic cohesion, our model-based clustering analysis of scores from the first five PC axes inferred 
only one inclusive cluster under the diagonal multivariate normal model, which had the lowest BIC value. 
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FIGURE 4. Scatterplots of factor scores from the first three principal component (PC) axes (A) and canonical scores from 
the first three canonical axes of a linear discriminant analysis (B) of multivariate phenotypic data from Crotalus scutulatus 
specimens. Symbol colors represent Mohave-Sonoran (green), Chihuahuan (yellow), Central Plateau (orange), and Crotalus 
scutulatus salvini (black) genetic lineages. 
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TABLE 1. Selected phenotypic characters for the genetic lineages of Crotalus scutulatus. For meristic characters, range is 
followed by mode, except for C. salvini due to small sample size or where there are multiple modal values, we presented 
the median (shown in bold). For snout vent length (SVL), range in mm is followed by the average. Values for males are 
shown above those for females. 


Characters 

C. scutulatus 

Mohave-Sonoran 

n = 167 

C. scutulatus 

Chihuahuan 

n= 114 

C. scutulatus 

Central Mexican 

Plateau 

n = 49 

C. scutulatus salvini 

n = 17 

Prefrontals 

4-20 (8) 

5-13 (6,8) 

6-12 (6) 

5-7(6) 


5-16(9,11) 

5-14 (8) 

5-11(6) 

6-7(6) 

Interrictals 

18-30 (26) 

18-28 (24) 

20-25 (22) 

18-23 (21) 


21-28 (25) 

20-26 (23) 

19-25 (24) 

19-21 (20) 

Supralabials 

12-17(15) 

13-17(15) 

12-17(14) 

12-15 (13.5) 


13-16(15) 

13-18(15) 

12-16(14) 

12-15 (13) 

Infralabials 

13-19(16) 

14-18 (16) 

14-17(15) 

14-16 (14.5) 


14-17(15) 

14-18(16) 

13-17(15) 

13-16 (14.5) 

Ventrals 

169-189 (176) 

164-181 (173) 

162-176 (165) 

164-180 (169) 


175-190 (184) 

167-183 (177) 

166-182 (168) 

168-182 (177) 

Subcaudals 

19-32 (25) 

16-29 (25) 

13-26 (20) 

16-28 (24) 


14-24(19) 

15-25(18) 

11-21 (17) 

16-21 (17) 

Rattle fringe 

10-13 (12) 

9-14(12) 

10-13 (12) 

10-14 (12) 


10-14 (11) 

9-13(12) 

10-14(12) 

11-14 (12) 

Dorsal scale rows 

23-28 (25) 

24-28 (25) 

25-28 (25) 

22-30 (25) 


22-28 (25) 

24-27 (25) 

24-29 (25) 

24-26 (25) 

Dorsal body blotches 27-43 (37) 

27^43 (36) 

25-39 (37) 

28-33 (30) 


32-42 (36) 

31-47 (36) 

30-42 (33) 

29-35 (33) 

Tail bands 

3-8 (5) 

2-6(4) 

2-6 (4) 

3-7(5) 


2-6 (4) 

2-6 (3) 

2-5 (2,3) 

3—4 (3) 

SVL 

478-1066(752) 

442-997 (722) 

370-929 (624) 

645-1067(761) 


478-926 (702) 

322-853 (616) 

442-808 (610) 

470-833 (757) 


Principal component axis 1 (42.1% of variance explained) reflected correlated variation in morphometric traits, 
for which all factor loadings weighted negatively, along with meristic traits SubC, Irict, and TB (Fig. 5; see Appen¬ 
dix 2 for character abbreviations). By contrast, PC 2 (11.5% of variance explained) reflected correlated variation 
in primarily meristic traits, with high values negatively correlated with this axis (Fig 5). Four characters, FIL, HW, 
MBL, and MBI were positively correlated with PC2 with moderately high factor loadings, and this axis also showed 
the greatest differences in mean values across genetic lineages as indicated by the highest /-value from ANOVA 
comparisons (Fig. 5; Table 2). The PC 3 (7.7% of variance explained) reflected aspects of sexual dimorphism that 
were consistent across genetic lineages, with CL, SubC, DBB, and TB positively correlated, and DSR, Slab, Ilab, 
and Irict, inversely correlated with this axis, respectively. By contrast, PC 4 (7.2% of variance explained) reflected 
sexually dimorphic traits that varied across genetic lineages (lineage and sex interaction), with BL, Vent, PreF, and 
DBB positively correlated, and CL, SubC, Slab, Ilab, and TB negatively correlated along this axis. Finally, variation 
in PC 5 (5% of variance explained) did not reflect a clear pattern with respect to phenotypic variation across lineages 
or between sexes, and we interpret this axis to represent primarily “noise” variation consistent across both lineages 
and sexes. 

Two-way ANOVAs with lineage and sex as factors inferred a general pattern of significant differences between 
lineage and sexes in mean PC factor scores, with minimal interaction effects between factors (Fig. 5; Table 2). Least 
squares linear regression analyses of the first five PC axes against latitude and longitude indicated a general pattern 
of slopes that differed significantly from zero (Table 3). With the exception of PC 2, latitude and longitude generally 
explained little variation in PC factor scores (i.e., adjusted /?--values < 0.2); however, PC 2 demonstrated a relatively 
strong relationship between factor scores and latitude ( R 2 = 0.43) as well as between factor scores and longitude (R 2 
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= 0.31; Fig. 6). Because PC 4 reflected variation among sexually dimorphic traits and lineages simultaneously, this 
axis also had relatively high /C-valucs for the effects of latitude and longitude (Table 3). Mantel tests also revealed 
an overall morphological pattern consistent with “isolation-by-distance” (R = 0.17, p < 0.0001). In our analysis of 
relative head length among lineages, we did not reject the null hypothesis of equal slopes (p = 0.21), and also found 
significant differences in relative head size (i.e., intercepts of FIL on BL) among genetic lineages of C. scutulatus, 
with specimens of C. s. salvini and the Central Plateau lineages showing greater relative head size than specimens 
from the northern two lineages (F= 2556.87, df= 1,3, p < 0.0001; Fig. 7). 
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FIGURE 5. Boxplots of factor scores from the first five principal component (PC) axes for each lineage separated by sex (left 
panel), and PC factor loadings for each character on respective PC axes (right panel). Lineage abbreviations are as follows: 
MOS (Mohave-Sonoran), CHI (Chihuahuan), CMP (Central Plateau), SAL (Crotalus scutulatus salvini). Character abbrevia¬ 
tions are provided in Appendix 2. 
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FIGURE 6. Principal component (PC) factor scores for axis 2 regressed on latitude and longitude, showing clinal variation of 
phenotypes within Crotalus scutulatus (left panel), and photographs in life of C. scutulatus individuals representative of each 
genetic lineage (right panel). Symbol colors represent Mohave-Sonoran (green), Chihuahuan (yellow), Central Plateau (orange), 
and C. s. salvini (black) lineages. Photo credits from top to bottom: JMM (south of Wellton, Yuma Co., Arizona), JMM (west of 
Davis Mountains, Jeff Davis Co., Texas), Michael Price (Huizache, San Luis Potosi, Mexico), and Carlos Hernandez-Jimenez 
(Jalapasco Tepayahualco, Puebla, Mexico). 

TABLE 2. Two-way ANOVA results with Lineage and Sex as factors and factor scores from the first five PC axes as the 


dependent variable. 

Dependent 

Lineage 



Sex 



Lineage:Sex 


Variable 

DF 

Rvalue 

/>-value 

DF 

Rvalue 

p-value 

DF 

Rvalue 

/»-value 

PC 1 

3 

10.922 

<0.0001 

1 

35.240 

<0.0001 

3 

3.616 

0.0139 

PC 2 

3 

57.060 

<0.0001 

1 

7.212 

0.0078 

3 

0.972 

0.4068 

PC 3 

3 

7.411 

<0.0001 

1 

25.362 

<0.0001 

3 

2.528 

0.0581 

PC 4 

3 

27.850 

<0.0001 

1 

244.493 

<0.0001 

3 

0.592 

0.621 

PC 5 

3 

1.821 

0.144 

1 

0.581 

0.447 

3 

0.654 

0.581 
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TABLE 3. Least squares linear regressions of factor scores from the first five principal component (PC) axes against 
latitude and longitude. 


Source 


B 

SEB 

t-value 

p-value 

Adjusted R 2 


Latitude 

-0.161 

0.046 

-3.499 

0.0006 

0.044 

PC 1 


Longitude 

0.088 

0.036 

2.419 

0.0163 

0.020 


Latitude 

-0.251 

0.019 

-13.44 

<0.0001 

0.426 

PC 2 


Longitude 

0.165 

0.016 

10.35 

<0.0001 

0.305 

PC 3 

Latitude 

0.048 

0.020 

2.423 

0.0161 

0.020 

Longitude 

-0.054 

0.015 

-3.498 

0.0006 

0.044 


Latitude 

0.091 

0.019 

4.911 

<0.0001 

0.087 

PC 4 


Longitude 

-0.097 

0.014 

-7.002 

<0.0001 

0.166 

PC 5 

Latitude 

0.027 

0.016 

1.676 

0.095 

0.007 

Longitude 

-0.019 

0.013 

-1.479 

0.140 

0.005 



FIGURE 7. Scaling relationship of head length (log mm) and body length (log mm) for each genetic lineage of Crotalus 
scutulatus. Lineage abbreviations are as follows: MOS (Mohave-Sonoran), CHI (Chihuahuan), CMP (Central Plateau), SAL 
(Crotalus scutulatus salvini). 


Discussion 

Our results demonstrate that Crotalus scutulatus exhibits variation in morphology that follows a cline along the 
major axis of its expansive geographic distribution. Relative to southeastern specimens, Mohave rattlesnakes from 
the northwestern part of their range had relatively high values for CL, Vent, SubC, DSR, Slab, Ilab, PreF, Irict, DBB, 
and TB, and relatively low values for F1L, FIW, MBL, and MBI (Figs. 5-6). Although this finding is consistent with 
early interpretations of variation in C. scutulatus (e.g., Klauber 1930; Gloyd 1940), no previous studies have evalu¬ 
ated phenotypic variation across the distribution of this taxon in an explicit statistical and geographical context. The 
major genetic lineages differed in mean values of multiple traits with respect to PC factor scores (Table 2). Overall 
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phenotypic variation, however, was continuous and lacked any discrete gaps or shifts in multivariate distributions 
across genetic lineage boundaries. Importantly, the results from model-based clustering indicate that all C. scutula¬ 
tus sampled are included in a single cohesive phenotypic cluster that was diagonal multivariate normal. 

Along with previously published genetic data, our phenotypic dataset provides additional insight into patterns 
of variation and divergence relevant to species delimitation in C. scutulatus. Hypothesis A posits that each of the 
four genetic lineages inferred by Schield et al. (2018) should be elevated to species (Fig. 1). This hypothesis ac¬ 
knowledges past population fragmentation and divergence within an ancestral lineage that temporally coincides 
with Pleistocene glacial cycles, and emphasizes the signature of these demographic events in genetic data (Schield 
et al. 2018). However, this hypothesis does not necessarily capture the complex demographic history and patterns of 
gene flow between genetic lineages, as demonstrated by evaluations of allele frequency spectra between pairwise C. 
scutulatus lineages, which indicate secondary contact or migration at all boundaries within major genetic lineages. 
In the present study, phenotypic data are consistent with expectations of past and/or ongoing gene flow between 
all adjacent C. scutulatus lineages, and in terms of morphological distinction or diagnosis, distinguishing species 
within a four-species model would be untenable. 

Hypothesis B, or the “expanded salvini ’ concept, posits that species boundaries coincide with the deepest ge¬ 
netic lineage divergence, whereby C. scutulatus would consist of Mohave-Sonoran and Chihuahuan populations, 
and C. salvini would comprise all populations from the Central Plateau southward (Fig. 1). One of the merits of 
this hypothesis is that Central Plateau and C. s. salvini populations together form a clade in phylogenetic infer¬ 
ences (Schield et al. 2018). However, Schield et al. (2018) also inferred relatively low levels of nuclear allelic dif¬ 
ferentiation between Central Plateau and Chihuahuan lineages, which could indicate relatively high levels of gene 
flow between these groups, or could be an artifact of the large geographic area without molecular samples between 
the Mohave-Sonoran and Chihuahuan groups. Morphological differentiation across the Central Plateau and Chi¬ 
huahuan clade boundary was also minimal, and the low average posterior probability of classification of Central 
Plateau specimens primarily was due to high misclassification of specimens between these lineages. Hypothesis C 
is similar to Hypothesis B in that it is consistent with a combined interpretation of mtDNA and nuclear coalescent 
species trees, but also suffers the same drawbacks to Hypothesis B in that it is not supported by phenotypic data and 
shows both low nuclear differentiation and a relatively strong signal for historical gene flow between Chihuahuan 
and Central Plateau populations. The “expanded salvini ” concepts are somewhat supported by the observation that 
Central Plateau and Crotalus s. salvini populations tend to have proportionately larger heads than do snakes from 
the northern lineages. However, relative HL contributes to the overall pattern of clinal variation in C. scutulatus 
phenotypes (Figs. 5-6), and shows high variation within and among populations (Fig. 7). 

Hypothesis D is based partially on the legacy taxonomy of the past century, when both subspecies of C. scu¬ 
tulatus were formalized based on putative phenotypic distinctions (Klauber 1930, 1972; Gloyd 1940; Fig. 1). Our 
results are in general agreement with phenotypic data presented by Gloyd (1940) and Klauber (1930, 1972), but 
these early assessments were based largely on comparisons of specimens from the extreme northern and southern 
parts of the distribution of C. scutulatus, with relatively few specimens available to confirm clinal variation from the 
large gap across central Mexico. Individuals from the extreme northwestern and southeastern parts of C. scutulatus 
distribution are considerably different from each other in overall coloration and pattern (Fig. 6), with C. s. salvini 
having border-free chocolate dorsal blotches against a buff to pink background color, and often with striking jade- 
colored irises. However, our more expansive dataset including numerous specimens from central Mexico confirms 
the overall pattern of clinal variation predicted by Gloyd (1940), and further diminishes the diagnosability of C. s. 
salvini as a species distinct from C. s. scutulatus. In addition to a gradual increase in values of most meristic traits 
from south to north, we also note that the condition of distinct pale blotch borders increases in frequency from south 
to north across the distribution of C. scutulatus, rather than being diagnostic of C. s. scutulatus (Fig. 3). Although C. 
s. salvini constitutes a distinct lineage based on mitochondrial gene tree and nuclear coalescent species tree analy¬ 
ses, considering C. s. salvini an independent lineage/species from the Central Plateau lineage had the lowest support 
of all four lineages tested based on nuclear data in Schield et al. (2018). Based on a large nuclear SNP dataset, the 
demographic history at the boundary between C. 5. scutulatus and C. 5. salvini was inferred as “divergence with 
ancestral symmetrical migration” (Schield et al. 2018), suggesting pervasive gene flow or very recent divergence. 
However, as noted by Schield et al. (2018), these mitochondrial/nuclear discrepancies do not preclude the possibil¬ 
ity that C. s. salvini could represent a peripheral isolate population in the early process of speciation that has yet to 
complete lineage sorting. Unfortunately, there is a dearth of material available from the southern Central Plateau 
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region (e.g., Queretaro, Hidalgo) despite recent field efforts, and it is unclear if this reflects sampling bias or an 
actual recent distributional hiatus between C. s. scutulatus and C. s. salvini. 

Collectively, phenotypic data in context with previously published molecular studies suggest that Hypothesis E, 
the single species model, is the most parsimonious delimitation of species within C. scutulatus (Fig. 1). As expected, 
phenotypic data were generally congruent with nuclear data and support the notion that isolated geographic lineages 
of C. scutulatus currently, or until recently, have experienced secondary contact and renewed gene flow. Further¬ 
more, our finding of continuous clinal variation within C. scutulatus emphasizes the importance of morphological 
data in robust species delimitation, as broadly expanded genetic datasets are likely to reveal complex demographic 
scenarios that can obscure a straightforward delineation of lineages. In such contexts, morphological data can help 
identify pragmatic taxonomic solutions where genetic data alone do not provide a clear solution. The juxtaposition 
of morphological and genetic datasets within C. scutulatus illustrates an interesting case of a wide-ranging species 
that harbors substantial genetic diversity, and that has experienced a unique demographic history involving exten¬ 
sive gene flow between cryptic genetic lineages, which we conclude has likely influenced the observed pattern of 
clinal variation in morphology. This perception of C. scutulatus conforms to the general lineage concept of species 
as metapopulation lineages (de Queiroz 1998, 2007), and exemplifies how both fine-scale molecular data and phe¬ 
notypic data can now delimit the boundaries of these (meta)populations within cohesive, widespread species. 
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APPENDIX 1. Specimens examined 


Crotalus scutulatus. MEXICO (« = 121): AGUASCALIENTES: CAS 87400; ENCB 13451; MVZ 275536; TCWC 38569; 
UTA 18360. CHIHUAHUA: BYU 13871-72, 15313-14, 15321, 15344, 15349-51, 15678, 17108, 17113, 19133, 21717; KU 
35093; 45339-45, 62865, 75643; MVZ 68913, 71036-37, 71040, 71050, 73115, 73117, 73123, 84509; TNHC 101199, 105694; 
UTA 4554, 12587, 17932, 58939. COAHUILA: CAS-SU 27236; FMNH 44317; KU 39952; TCWC 38991, 54241; TNHC 
29870, 33887, 112044. DURANGO: BYU 40051-52; CAS 120880; FMNH 103000; ENCB 16074; KU 39955; MVZ 143532; 
MZFC 17996; TCWC 43327; TNHC 29761. GUANAJUATO: MVZ 275539. JALISCO: KU 95973-75, 102987, 102990-92. 
NUEVO LEON: KU 187318; SDNHM 57004; TCWC 36012, 54242, 60785-86; TNHC 29866, 29871, 112041; UTA 4595. 
PUEBLA: CNAR 6970; FMNH 106061; ENCB 1116; MVZ 275548, 275552 (MZFC 26289); MZFC 4305, 25208. SAN LUIS 
POTOSI: CAS 169831; ENCB 108, 8151; MVZ 275546^17. SONORA: CNAR 48, 6978; MVZ 81284-85; SDNHM 2327, 
37635, 49711. TAMAULIPAS: SDNHM 6573. TLAXCALA: CNAR 6218. VERACRUZ: BYU 42717, 42726; KU 23750, 
26688; MVZ 275543-45; MZFC 2008, 25166. ZACATECAS: CAS 96077; KU 29489, 29490, 39956, 45000; MVZ 143533, 
275549-50; SDNHM 49713; TCWC 56448; UTA 2715. 

UNITED STATES (n = 226): ARIZONA: Cochise Co.: BYU 5406; FMNH 207793; MVZ 209102, 209105, 215712, 226243, 
229828, 229829, 229831, 229836; TCWC 85929; TNHC 60138, 105691, 105692, 112040; UTA 10024, 10025, 35042, 38798, 
38800. Coconino Co.: SDNHM 2105-08, 2130-32. Graham Co.: CAS 170506; SDNHM 40618. La Paz Co.: MVZ 233686, 
233690-93. Maricopa Co.: CAS 65089, 65699; MVZ 129357, 129359, 233684; SDNHM 25556-57, TNHC 58831. Mohave 
Co.: MVZ 275561-62, 277717; SDNHM 867, 1145, 1820, 3484, 17305,29702,31779,39023,41584. Pima Co.: CAS 192717; 
FMNH 207794; MVZ 196873, 233673; SDNHM 2766, 23235; UTA 14465, 32336-37, 32343. Pinal Co.: CAS 192742; SD¬ 
NHM 2767, 17074. Santa Cruz Co.: CAS 80762; MVZ 76672; SDNHM 17948, 20926. Yavapai Co.: CAS 63881-83, 65111, 
65693, 156178; MVZ 70299, 97145^16, 233696-97; SDNHM 2025, 2044-46, 2080, 2122-24, 2182, 3039, 3463, 5449, 5566- 
67, 9302, 17628-30. YUMA Co.: SDNHM 2684, 2708; UTA 45169. CALIFORNIA: Kern Co.: CAS 91629; KU 62904; MVZ 
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31734, 137600, 193456-58, 193460, 233699. Los Angeles Co.: BYU 42690, 42708; MVZ 11430, 275565; SDNHM 838, 
865, 1656-57, 2409, 2810, 20118, 33688, 37967, 43609. San Bernardino Co.: BMNH 1987.1231; BYU 41724; CAS 228094, 
250115, 259911; KU 31352; MVZ 18034, 26666, 176158, 193464, 233782-83; SDNHM 2605, 10549. NEVADA: Clark Co.: 
BYU 52748, 62939, 63056; CAS 156168; MVZ 275567, 277732, 277738-39; SDNHM 4396. Lincoln Co.: MVZ 14304. NEW 
MEXICO: Hidalgo Co.: TCWC 56318; TNHC 104537; UTA51371, 53919-20, 53924-25. TEXAS: Brewster Co.: KU 13782, 
61326, 177079-82, 319075, 336523, 336525, 341047; TCWC 4748, 16148,16150,27534-35,29404, 33377, 36356-57, 36906, 
40098, 40100, 64342, 64556, 86017, 88019, 103248; TNHC 33886, 85217, 89444, 92383; UTA 53729. Culberson Co.: TCWC 
33376, 77638; UTA 44383. Hudspeth Co.: TCWC 18394; TNHC 87883. Jeff Davis Co.: CAS 170519; TCWC 72829, 84534, 
91848; TNHC 66531, 66539, 66881-82, 87880, 89782, 89786, 89855, 89857, 97400-01; UTA 3901, 53945, 61674. Pecos 
Co.: SDNHM 916; TNHC 89831; UTA 44376. Presidio Co.: KU 336524; TCWC 27533, 80604, 90727; TNHC 66436, 66530, 
66532-33. UTAH: Washington Co.: BYU 35963, 42699, 42849. 


APPENDIX 2. Morphological characters 

1. Body length, measured in mm (BL); 

2. Tail length, measured in mm (CL); 

3. Head length measured to 0.1 mm from face of the rostral to the angle of the quadrate and mandible (HL); 

4. Head width measured to 0.1 mm at widest point (HW); 

5. Snout length measured to 0.1 mm from anterior edge of the eye to the center of the rostral plate (SL); 

6. Eye-to-pit distance measured to 0.1 mm from the shortest distance between anterior edge of eye to posterior edge of pit 
organ (EPD); 

7. Intersupraocular distance measured to 0.1 mm from the lateral edge of each supraocular scale (ISD); 

8. Eye diameter measured to 0.1 mm along the midline of eye from anterior to posterior edge (Eye); 

9. Width of the basal rattle segment, measured in 0.1 mm along dorsal-ventral axis of rattle (BRW); 

10. Middle blotch length measured in 0.1 mm from the anterior to posterior border of the middle dorsal blotch along ver¬ 
tebral line (MBL); 

11. Posterior middle blotch interspace measured in 0.1 mm as the length of the interspace immediately posterior to the 
middle dorsal blotch along vertebral line (MBI); 

12. Number of ventral scales (Vent); 

13. Number of subcaudal scales (SubC); 

14. Number of dorsal scale rows at midbody (DSR); 

15. Number of rattle-fringe scales (RFS); 

16. Average number of supralabial scales on the left and right side of the head (SLab); 

17. Average number of infralabial scales on the left and right side of the head (ILab); 

18. Number of scales in the prefrontal region, including all scales anterior to the supraoculars but excluding the intemasals 
(PreF); 

19. Number of interrictal scales (Irict)—counted as the number of dorsal head scale rows between each rictus, excluding 
ultimate supralabials; 

20. Number of dorsal body blotches (DBB)—interconnected blotches counted as one except for three or more, which were 
counted separately; 

21. Number of tail bands (TB)—interconnected bands counted as one except for three or more, which were counted sepa¬ 
rately; 

22. Condition of dorsal blotch borders (SBB)—described as distinct, indistinct, or absent. 
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